#define MECPlotting_cxx
#include "MECPlotting.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <string>
#include <list>
#include <functional>
#include <numeric>
#include <TClonesArray.h>
#include <TRefArray.h>
#include <TMath.h>

#include <TSystem.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TF1.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TTree.h>
#include "TDirectory.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TString.h"
#include <TLegend.h>
#include "TLine.h"


#define TRUE 0
#define FALSE 1
#define PROTON 2212
#define NEUTRON 2112
#define MUON 13


using namespace std;
using std::cout;
using std::endl;

TList* addHistogram1D(TList*, const char*, const char*, int, double, double);
TList* addHistogram2D(TList*, const char*, const char*, int, double, double, int, double, double);
TList* initialiseHists(TList*);

//void MECPlotting::Loop()
int analysis(int NEvents=0, TString fileName="prod6_analysis.list", TString outSaveName="histSave.root"){ 
{
  //Specify Error Flag
  int ISERROR = FALSE;
  
  //Initialise Histogram list
  TList* histList = new TList();

  //create hist templates  
  histList = initialiseHists(histList);  
  
  // Get list of files to run over. 
  std::ifstream inputFile(fileName.Data(), ios::in);
  
  // Declare a TChain for the data we care about
  TChain *dataChain = new TChain("tree");
  
  // Check if the file exists.
  if (!inputFile.is_open())
  {
    std::cout << "ERROR: File list not found!" << std::endl;
    std::cout << " - File should contain list of files to be processed." << std::endl;
    return(1);
  }
  else
  {
    std::string curFileName;
    // Add the input files to the TChains.
    while(getline(inputFile,curFileName))
    {
      dataChain->Add(curFileName.c_str());
    }
  }

  std::cout << "Got input file(s)." <<std:: endl;
  
  ISERROR = fillHists(NEvents, dataChain, histList);
  
  //Wrap Up:

  TFile histSave(outSaveName,"RECREATE");

  histList->Write();
  
  //Canvas creation and combined histograms goes here
  
  
    histSave.Close();

    std::cout << "***** Dubug info ***** " << std::endl;
    std::cout << "GLOB_i is " <<  GLOB_i << std::endl;
    std::cout << "GLOB_VTX is " <<  GLOB_VTX << std::endl;
    std::cout << "GLOB_FGD is " <<  GLOB_FGD << std::endl;
    std::cout << "GLOB_MEC is " <<  GLOB_MEC << std::endl;
    std::cout << "GLOB_MEC_FGD is " <<  GLOB_MEC_FGD << std::endl;
    std::cout << "GLOB_NOTMEC is " <<  GLOB_NOTMEC << std::endl;
    std::cout << "GLOB_NOTMEC_FGD is " <<  GLOB_NOTMEC_FGD << std::endl;

    if (ISERROR == TRUE)
    {
      cout << "Error: Exiting now ..." << endl;
      return(1);
    }
    else
    {
      cout << "Successful Run, Exiting now ..." << endl;
      return(0);
    }
  
}


TList* addHistogram1D(TList* histList, const char* name, const char* title, int nbins, double min,double max){
  TH1* hist = new TH1D(name, title, nbins, min, max);
  histList->Add(hist);
  return(histList);
}



TList* addHistogram2D(TList* histList, const char* name, const char* title, int nbins_x, double min_x, double max_x, int nbins_y, double min_y, double max_y){
  TH2* hist = new TH2D(name, title, nbins_x, min_x, max_x, nbins_y, min_y, max_y);
  histList->Add(hist);
  return(histList);
}

TList* initialiseHists(TList* histLis){
  histList = addHistogram1D(histList, "EnuRecHistFromP","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14000);
  histList = addHistogram1D(histList, "EnuRecHistFromP_t","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_MEC_t","Histogram of Recontructed MEC Target Mass inside FGD",600,-2000,6000);
  histList = addHistogram1D(histList, "EnuRecHistFromP_NotMEC","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14000);
  histList = addHistogram1D(histList, "EnuRecHistFromP_NotMEC_t","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_NotMEC_t","Histogram of Recontructed not MEC Target Mass inside FGD",600,-2000,6000);

  histList = addHistogram1D(histList, "ReconTargetMassFromP_MEC","Histogram of Recontructed MEC Target Mass inside FGD",300,-500,3000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_CCQE","Histogram of Recontructed CCQE Target Mass inside FGD",300,-500,3000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_RES","Histogram of Recontructed RES Target Mass inside FGD",300,-500,3000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_OTHER","Histogram of Recontructed OTHER Target Mass inside FGD",300,-500,3000);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_NotMEC","Histogram of Recontructed not MEC Target Mass inside FGD",300,-500,3000);



  histList = addHistogram1D(histList, "EnuRecHist_comb","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14000);
  histList = addHistogram1D(histList, "ReconTargetMass_comb","Histogram of Recontructed MEC Target Mass inside FGD",600,-2000,6000);


  histList = addHistogram1D(histList, "nOfFinalStateNucleons_MEC","Number of Nucleons in final state of MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfFinalStateProtons_MEC","Number of Protons in final state of MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfFinalStateNeutrons_MEC","Number of Neutrons in final state of MEC event",10,0,10);

  histList = addHistogram1D(histList, "nPIDS_MEC","Number of recon PIDS in MEC Event",10,0,10);
  histList = addHistogram1D(histList, "nOfReconProtons_MEC","Number of recon Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfReconMuons_MEC","Number of recon Muons in MEC event",10,0,10);

  histList = addHistogram1D(histList, "nOfTruthProtons_MEC","Number of true Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfTruthMuons_MEC","Number of true Muons in MEC event",10,0,10);

  histList = addHistogram1D(histList, "nOfTruthProtonsTPC_MEC","Number of true Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfTruthMuonsTPC_MEC","Number of true Muons in MEC event",10,0,10);

  histList = addHistogram2D(histList, "nOfTrueProtonsVsReconProtons","nOfTrueProtonsVsReconProtons",10,0,10,10,0,10);
  histList = addHistogram2D(histList, "nOfTrueProtonsTPCVsReconProtons","nOfTrueProtonsVsReconProtons",10,0,10,10,0,10);
  histList = addHistogram2D(histList, "nOfTrueProtonsVsProtonMom","nOfTrueProtonsVsReconProtonMom",10,0,10,50,0,12000);

  histList = addHistogram2D(histList, "EnuRecVsReconTargetMass_MEC","Histogram of Recontructed Neutrino Energy vs Recon Target Mass (MEC) ",300,0,14000,600,-4000,4000);
  histList = addHistogram2D(histList, "EnuRecVsReconTargetMass_NotMEC","Histogram of Recontructed Neutrino Energy vs Recon Target Mass (Not MEC) ",300,0,14000,600,-4000,4000);

  histList = addHistogram1D(histList, "nuAngle","Histogram of incoming neutrino angle",100,0,3);
  histList = addHistogram1D(histList, "nuAngle_t","Histogram of incoming neutrino angle",100,0,3);


  histList = addHistogram1D(histList, "pP_tHistInFGD","Histogram of MEC Proton Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "pP_tHistInFGD_NotMEC","Histogram of Not MEC Proton Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "mP_tHistInFGD","Histogram of MEC muon Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "mP_tHistInFGD_NotMEC","Histogram of Not MEC muon Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "difP_tHistInFGD","Histogram of MEC dif between muon and proton Pt",50,-2000,2000);
  histList = addHistogram1D(histList, "difP_tHistInFGD_NotMEC","Histogram of Not MEC dif between muon and proton Pt",50,-2000,2000);

  histList = addHistogram1D(histList, "pP_tHistInFGD_t","Histogram of MEC Proton Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "pP_tHistInFGD_NotMEC_t","Histogram of Not MEC Proton Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "mP_tHistInFGD_t","Histogram of MEC muon Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "mP_tHistInFGD_NotMEC_t","Histogram of Not MEC muon Transverse Momentum (final final)",50,-1000,3000);
  histList = addHistogram1D(histList, "difP_tHistInFGD_t","Histogram of MEC dif between muon and proton Pt",50,-2000,2000);
  histList = addHistogram1D(histList, "difP_tHistInFGD_NotMEC_t","Histogram of Not MEC dif between muon and proton Pt",50,-2000,2000);

  histList = addHistogram1D(histList, "deltaPhiPt","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_CCQE","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_RES","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_OTHER","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_NotMEC","Histogram of deltaPhiPt",100,0,3);

  histList = addHistogram1D(histList, "deltaPhiPt_t","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_NotMEC_t","Histogram of deltaPhiPt",100,0,3);

  histList = addHistogram1D(histList, "pPHistAfterCuts_MEC","Histogram of  MEC Proton Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistAfterCuts_CCQE","Histogram of  CCQE Proton Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistAfterCuts_RES","Histogram of  CCRES Proton Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistAfterCuts_OTHER","Histogram of  OTHER Proton Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistAfterCuts_NotMEC","Histogram of  Not MEC Proton Momentum after cuts",70,0,3000);

  histList = addHistogram1D(histList, "pPHistBeforeCuts_MEC","Histogram of  MEC Proton Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_CCQE","Histogram of  CCQE Proton Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_RES","Histogram of  CCRES Proton Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_OTHER","Histogram of  OTHER Proton Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_NotMEC","Histogram of  Not MEC Proton Momentum before cuts",70,0,3000);

  histList = addHistogram1D(histList, "mPHistAfterCuts_MEC","Histogram of  MEC Muon Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistAfterCuts_CCQE","Histogram of  CCQE Muon Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistAfterCuts_RES","Histogram of  CCRES Muon Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistAfterCuts_OTHER","Histogram of  OTHER Muon Momentum after cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistAfterCuts_NotMEC","Histogram of  Not MEC Muon Momentum after cuts",70,0,3000);

  histList = addHistogram1D(histList, "mPHistBeforeCuts_MEC","Histogram of  MEC Muon Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_CCQE","Histogram of  CCQE Muon Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_RES","Histogram of  CCRES Muon Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_OTHER","Histogram of  OTHER Muon Momentum before cuts",70,0,3000);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_NotMEC","Histogram of  Not MEC Muon Momentum before cuts",70,0,3000);

  return(histList);
 
}
